#include "Solver.h"


//====================================
// AMG Solver
//====================================

void AMG(HYPRE_ParCSRMatrix A0, HYPRE_ParVector b0, HYPRE_ParVector x0){
    
    if(myrank == 0) printf("\nUsing AMG solver...\n");

    HYPRE_Solver solver;
    HYPRE_BoomerAMGCreate(&solver);
    HYPRE_BoomerAMGSetTol(solver, 1.0e-6);
    HYPRE_BoomerAMGSetMaxIter(solver, 2000);
    HYPRE_BoomerAMGSetup(solver, A0, b0, x0);
    HYPRE_BoomerAMGSolve(solver, A0, b0, x0);

    int Iteration;
    double res;

    HYPRE_BoomerAMGGetNumIterations(solver, &Iteration);
    HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &res);

    if(myrank == 0){
        printf("\n");
        printf("\nIterations = %d\n", Iteration);
        printf("\nFinal Relative Residual Norm = %e\n", res);
        }
    HYPRE_BoomerAMGDestroy(solver);
    }


//====================================
// CG Solver
//====================================

void CG(HYPRE_ParCSRMatrix A0, HYPRE_ParVector b0, HYPRE_ParVector x0){
    
    if(myrank == 0) printf("\nUsing CG solver...\n");

    HYPRE_Solver solver;
    HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
    HYPRE_PCGSetMaxIter(solver, 2000);
    HYPRE_PCGSetTol(solver, 1e-6);
    HYPRE_PCGSetTwoNorm(solver, 1);

    HYPRE_ParCSRPCGSetup(solver, A0, b0, x0);
    HYPRE_ParCSRPCGSolve(solver, A0, b0, x0);
	
	//printf("\ndone creating solver");

    int Iteration;
    double res;

    HYPRE_PCGGetNumIterations(solver, &Iteration);
    HYPRE_PCGGetFinalRelativeResidualNorm(solver, &res);

    if(myrank == 0){
        printf("\n");
        printf("\nIterations = %d\n", Iteration);
        printf("\nFinal Relative Residual Norm = %e\n", res);
        }
    HYPRE_ParCSRPCGDestroy(solver);
    }
    
    
//====================================
// PCG Solver with AMG Preconditioner
//====================================

void CG_AMG(HYPRE_ParCSRMatrix A0, HYPRE_ParVector b0, HYPRE_ParVector x0){
    
    if(myrank == 0) printf("\nUsing CG solver with AMG precondi...\n");

    HYPRE_Solver solver, precond;

    HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
    
    HYPRE_PCGSetTol(solver, 1.0e-6);
    HYPRE_PCGSetMaxIter(solver, 2000);
    HYPRE_PCGSetTwoNorm(solver, 1);

    HYPRE_BoomerAMGCreate(&precond);
    HYPRE_BoomerAMGSetCoarsenType(precond, 6);
    HYPRE_BoomerAMGSetRelaxType(precond, 6);
    HYPRE_BoomerAMGSetNumSweeps(precond, 1);

    HYPRE_BoomerAMGSetTol(precond, 0.0);
    HYPRE_BoomerAMGSetMaxIter(precond, 1);
    
    HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup, precond);
    
    HYPRE_ParCSRPCGSetup(solver, A0, b0, x0);
    HYPRE_ParCSRPCGSolve(solver, A0, b0, x0);
    
    int Iteration;
    double res;

    HYPRE_PCGGetNumIterations(solver, &Iteration);
    HYPRE_PCGGetFinalRelativeResidualNorm(solver, &res);

    if(myrank == 0){
        printf("\n");
        printf("\nIterations = %d\n", Iteration);
        printf("\nFinal Relative Residual Norm = %e\n", res);
        }
    HYPRE_ParCSRPCGDestroy(solver);
    HYPRE_BoomerAMGDestroy(precond);
    }
    
    
//====================================
// GMRES Solver with AMG Preconditioner
//====================================

int hypre_FlexGMRESModifyPCAMGExample(void *precond_data, int iterations, double rel_residual_norm){
    if(rel_residual_norm > .1){
        HYPRE_BoomerAMGSetNumSweeps(precond_data, 10);
        }
    else{
        HYPRE_BoomerAMGSetNumSweeps(precond_data, 1);
        }
    return 0;
    }    

void GMRES_AMG(HYPRE_ParCSRMatrix A0, HYPRE_ParVector b0, HYPRE_ParVector x0){
    
    if(myrank == 0) printf("\nUsing GMRES solver with AMG precondi...\n");

    int restart = 30;
    int modify  =  1;

    HYPRE_Solver solver, precond;

    HYPRE_ParCSRFlexGMRESCreate(MPI_COMM_WORLD, &solver);
    
    HYPRE_FlexGMRESSetTol(solver, 1.0e-6);
    HYPRE_FlexGMRESSetMaxIter(solver, 2000);
    HYPRE_FlexGMRESSetKDim(solver, restart);

    HYPRE_BoomerAMGCreate(&precond);
    HYPRE_BoomerAMGSetCoarsenType(precond, 6);
    HYPRE_BoomerAMGSetRelaxType(precond, 6);
    HYPRE_BoomerAMGSetNumSweeps(precond, 1);

    HYPRE_BoomerAMGSetTol(precond, 0.0);
    HYPRE_BoomerAMGSetMaxIter(precond, 1);
    
    HYPRE_FlexGMRESSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSolve, (HYPRE_PtrToSolverFcn) HYPRE_BoomerAMGSetup, precond);
    
    if (modify) 
        HYPRE_FlexGMRESSetModifyPC(solver, (HYPRE_PtrToModifyPCFcn) hypre_FlexGMRESModifyPCAMGExample);

    HYPRE_ParCSRFlexGMRESSetup(solver, A0, b0, x0);
    HYPRE_ParCSRFlexGMRESSolve(solver, A0, b0, x0);
    
    int Iteration;
    double res;

    HYPRE_FlexGMRESGetNumIterations(solver, &Iteration);
    HYPRE_FlexGMRESGetFinalRelativeResidualNorm(solver, &res);

    if(myrank == 0){
        printf("\n");
        printf("\nIterations = %d\n", Iteration);
        printf("\nFinal Relative Residual Norm = %e\n", res);
        }
    HYPRE_ParCSRFlexGMRESDestroy(solver);
    HYPRE_BoomerAMGDestroy(precond);
    }
    
    
//====================================
// PCG Solver with ParaSail Preconditioner
//====================================

void CG_ParaSails(HYPRE_ParCSRMatrix A0, HYPRE_ParVector b0, HYPRE_ParVector x0){
    
    if(myrank == 0) printf("\nUsing CG solver with ParaSails precondi...\n");
    
    int    sai_max_levels = 1;
    double sai_threshold  = 0.1;
    double sai_filter     = 0.05;
    int    sai_sym        = 1;
    int    restart        = 30;

    HYPRE_Solver solver, precond;

    HYPRE_ParCSRPCGCreate(MPI_COMM_WORLD, &solver);
    
    HYPRE_PCGSetTol(solver, 1.0e-6);
    HYPRE_PCGSetMaxIter(solver, 2000);
    HYPRE_PCGSetTwoNorm(solver, 1);

    HYPRE_ParaSailsCreate(MPI_COMM_WORLD, &precond);
    HYPRE_ParaSailsSetParams(precond, sai_threshold, sai_max_levels);
    HYPRE_ParaSailsSetFilter(precond, sai_filter);
    HYPRE_ParaSailsSetSym(precond, sai_sym);
    
    HYPRE_PCGSetPrecond(solver, (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSolve, (HYPRE_PtrToSolverFcn) HYPRE_ParaSailsSetup, precond);
    
    HYPRE_ParCSRPCGSetup(solver, A0, b0, x0);
    HYPRE_ParCSRPCGSolve(solver, A0, b0, x0);
    
    int Iteration;
    double res;

    HYPRE_PCGGetNumIterations(solver, &Iteration);
    HYPRE_PCGGetFinalRelativeResidualNorm(solver, &res);

    if(myrank == 0){
        printf("\n");
        printf("\nIterations = %d\n", Iteration);
        printf("\nFinal Relative Residual Norm = %e\n", res);
        }
    HYPRE_ParCSRPCGDestroy(solver);
    HYPRE_ParaSailsDestroy(precond);
    }
